!> author: 左志华
!> date: 2022-07-16
!>
!> Save particle data to file <br>
!> 存储粒子域信息
#ifdef REAL32
#define h5pt_writefileattrib_r8 h5pt_writefileattrib_r4
#define h5pt_writedata_r8 h5pt_writedata_r4
#endif
module sph_save_pif

    use h5part_m
    use seakeeping_time, only: now
    use sph_kinds, only: rk
    use sph_region, only: region_t
    use sph_terminal, only: term
    use sph_env, only: nums, hsml, dim

    integer(8), save :: nstep = 1_8 !! number of steps
                                    !! 步数
    logical :: save_acc     !! Save acceleration on/off <br>
                            !! 保存加速度

contains

    !> PIF handle <br>
    !> PIF 文件句柄
    subroutine handle_PIF(fid, proj_name, region, h5part_file, close)
        integer(8), intent(inout) :: fid        !! file id
                                                !! 文件句柄
        character(*), intent(in) :: proj_name   !! project name
                                                !! 项目名称
        type(region_t), intent(in) :: region
        character(len=*), intent(in), optional :: h5part_file
        logical, intent(in), optional :: close
        integer(8) :: status

        if (present(h5part_file)) then
            fid = h5pt_openw(h5part_file)
            if (fid == 0_8) then
                call term%error('file may be occupied by another program: '//h5part_file)
            end if

            status = h5pt_writefileattrib_string(fid, "CreateDate", now())
            status = h5pt_writefileattrib_string(fid, "Creator", "engine")
            status = h5pt_writefileattrib_string(fid, "Project", proj_name)
            status = h5pt_writefileattrib_i4(fid, "Nreal, Nvirt, Ntotal, Dimesion", &
                                             [nums(1:3), dim], 4_8)
            status = h5pt_writefileattrib_r8(fid, "HSML", [hsml], 1_8)

            status = h5pt_setnpoints(fid, int(nums(4), 8))
        else if (present(close)) then
            status = h5pt_close(fid)
        end if

    end subroutine handle_PIF

    !> Save particle data to file <br>
    !> 存储粒子域信息 (HDF5::H5part)
    subroutine save_PIF(fid, region, dim, n, t, istep)
        integer(8), intent(in) :: fid  !! H5part 文件单元号
        type(region_t), intent(in) :: region  !! 粒子宏观变量域
        integer, intent(in) :: dim, n  !! 粒子坐标维度及粒子数量
                                       !! @note 用以声明局部数组维度
        real(rk), intent(in) :: t  !! 时间
        integer(8), intent(inout) :: istep  !! 步数
        !@todo: 完成存储粒子域信息的功能

        integer(8) :: status; integer :: i
        real(rk) :: z(n, dim)

        status = h5pt_setstep(fid, istep)
        status = h5pt_writestepattrib(fid, "RealTime", [t], 1_8)

        z(:, 1) = 0.0_rk
        status = h5pt_writedata_r8(fid, "z", z(:, 1))  ! paraview 需要 z 变量

        ! @note 值复制输出，可能会存在性能问题，这是因为 h5part 以向量存储数据，方便可视化；Ti-SPH 以矩阵存储数据，
        ! 方便高效运算；因此需要将矩阵转换为向量。
        do concurrent(i=1:n)
            z(i, :) = region%particles%loc(:, i)
        end do

        status = h5pt_writedata_r8(fid, "x", z(:, 1))
        status = h5pt_writedata_r8(fid, "y", z(:, 2))

        do concurrent(i=1:n)
            z(i, :) = region%particles%vel(:, i)
        end do
        status = h5pt_writedata_r8(fid, "Vx", z(:, 1))
        status = h5pt_writedata_r8(fid, "Vy", z(:, 2))

        if (save_acc) then
            do concurrent(i=1:n)
                z(i, :) = region%particles%acc(:, i)
            end do
            status = h5pt_writedata_r8(fid, "Fx", z(:, 1))
            status = h5pt_writedata_r8(fid, "Fy", z(:, 2))
        end if

        status = h5pt_writedata_r8(fid, "Mass", region%particles%mass)
        status = h5pt_writedata_r8(fid, "Density", region%particles%rho)
        status = h5pt_writedata_r8(fid, "P", region%particles%p)
        status = h5pt_writedata_i4(fid, "Itype", region%particles%itype)
        istep = istep + 1_8

    end subroutine save_PIF

end module sph_save_pif
